#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX2
#include<string.h>
#include<iostream>
#include<fstream>
using namespace std;
const int N = 4100;
float nums[N][N];
int n;

long long head, tail, freq;
struct timespec sts, ets;
time_t dsec;
long dnsec;

ifstream fin("input.txt");
ofstream fout("output.txt");

void print_result(string str) {
	dsec = ets.tv_sec - sts.tv_sec;
	dnsec = ets.tv_nsec - sts.tv_nsec;
	if (dnsec < 0) {
		dsec--;
		dnsec += 1000000000ll;
	}
	ofstream of_result("/home/u149516/result.txt", ios::app);
	of_result << str << " with n= " << n << " time= " << dsec << "." << dnsec << "s" << endl;
	of_result.close();

}

void avx512_Gauss(int n) {
	for (int k = 1; k <= n; ++k) {
		__m512 vt = _mm512_set1_ps(nums[k][k]);
		int j = k + 1;
		for (; j + 16 <= n; j += 16) {
			__m512 va = _mm512_loadu_ps(&nums[k][j]);
			va = _mm512_div_ps(va, vt);
			_mm512_storeu_ps(&nums[k][j], va);
		}
		for (; j <= n; ++j)
			nums[k][j] /= nums[k][k];
		nums[k][k] = 1.0;
		for (int i = k + 1; i <= n; ++i) {
			int j = k + 1;
			__m512 vik = _mm512_set1_ps(nums[i][k]);
			for (; j + 16 <= n; j += 16) {
				__m512 vkj = _mm512_loadu_ps(&nums[k][j]);
				__m512 vij = _mm512_loadu_ps(&nums[i][j]);
				__m512 vx = _mm512_mul_ps(vik, vkj);
				vij = _mm512_sub_ps(vij, vx);
				_mm512_storeu_ps(&nums[i][j], vij);
			}
			for (; j <= n; ++j)
				nums[i][j] -= nums[i][k] * nums[k][j];
			nums[i][k] = 0;
		}
	}
}

void Gauss(int n) {
	for (int k = 1; k <= n; ++k) {
		for (int j = k + 1; j <= n; ++j)nums[k][j] /= nums[k][k];
		nums[k][k] = 1;
		for (int i = k + 1; i <= n; ++i) {
			for (int j = k + 1; j <= n; ++j)
				nums[i][j] -= nums[i][k] * nums[k][j];
			nums[i][k] = 0;
		}
	}
}

void sse128_Gauss(int n) {
	for (int k = 1; k <= n; ++k) {
		__m128 vt = _mm_set1_ps(nums[k][k]);
		int j = k + 1;
		for (; j + 4 <= n; j += 4) {
			__m128 va = _mm_loadu_ps(&nums[k][j]);
			va = _mm_div_ps(va, vt);
			_mm_storeu_ps(&nums[k][j], va);
		}
		for (; j <= n; ++j)
			nums[k][j] /= nums[k][k];
		nums[k][k] = 1.0;
		for (int i = k + 1; i <= n; ++i) {
			int j = k + 1;
			__m128 vik = _mm_set1_ps(nums[i][k]);
			for (; j + 4 <= n; j += 4) {
				__m128 vkj = _mm_loadu_ps(&nums[k][j]);
				__m128 vij = _mm_loadu_ps(&nums[i][j]);
				__m128 vx = _mm_mul_ps(vik, vkj);
				vij = _mm_sub_ps(vij, vx);
				_mm_storeu_ps(&nums[i][j], vij);
			}
			for (; j <= n; ++j)
				nums[i][j] -= nums[i][k] * nums[k][j];
			nums[i][k] = 0;
		}
	}

}

void avx256_Gauss(int n) {
	for (int k = 1; k <= n; ++k) {
		__m256 vt = _mm256_set1_ps(nums[k][k]);
		int j = k + 1;
		for (; j + 8 <= n; j += 8) {
			__m256 va = _mm256_loadu_ps(&nums[k][j]);
			va = _mm256_div_ps(va, vt);
			_mm256_storeu_ps(&nums[k][j], va);
		}
		for (; j <= n; ++j)
			nums[k][j] /= nums[k][k];
		nums[k][k] = 1.0;
		for (int i = k + 1; i <= n; ++i) {
			int j = k + 1;
			__m256 vik = _mm256_set1_ps(nums[i][k]);
			for (; j + 8 <= n; j += 8) {
				__m256 vkj = _mm256_loadu_ps(&nums[k][j]);
				__m256 vij = _mm256_loadu_ps(&nums[i][j]);
				__m256 vx = _mm256_mul_ps(vik, vkj);
				vij = _mm256_sub_ps(vij, vx);
				_mm256_storeu_ps(&nums[i][j], vij);
			}
			for (; j <= n; ++j)
				nums[i][j] -= nums[i][k] * nums[k][j];
			nums[i][k] = 0;
		}
	}

}
int main() {
	fin.open("input.txt");
	fin.close();
	fin.open("input.txt");
	fout.open("output.txt");
	fin >> n;
	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n; ++j)
			fin >> nums[i][j];
	fin.close();
	//timer
	timespec_get(&sts, TIME_UTC);
	Gauss(n);
	timespec_get(&ets, TIME_UTC);
	print_result("common x86 gauss");

		for (int i = 1; i <= n; ++i) {
			for (int j = 1; j <= n; ++j)
				fout << nums[i][j] << " ";
			fout << endl;
		}
	fout.close();
	memset(nums, 0, sizeof nums);

	fin.open("input.txt");
	fout.open("output.txt");
	fin >> n;
	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n; ++j)
			fin >> nums[i][j];
	fin.close();
	//timer
	timespec_get(&sts, TIME_UTC);
	sse128_Gauss(n);
	timespec_get(&ets, TIME_UTC);
	print_result("ss128 gauss");

		for (int i = 1; i <= n; ++i) {
			for (int j = 1; j <= n; ++j)
				fout << nums[i][j] << " ";
			fout << endl;
		}
	fout.close();
	memset(nums, 0, sizeof nums);

	fin.open("input.txt");
	fout.open("output.txt");
	fin >> n;
	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n; ++j)
			fin >> nums[i][j];
	fin.close();
	//timer
	timespec_get(&sts, TIME_UTC);
	avx256_Gauss(n);
	timespec_get(&ets, TIME_UTC);
	print_result("avx256 gauss");

	for (int i = 1; i <= n; ++i) {
		for (int j = 1; j <= n; ++j)
			fout << nums[i][j] << " ";
		fout << endl;
	}
	fout.close();
	memset(nums, 0, sizeof nums);

	
	fin.open("input.txt");
	fout.open("output.txt");
	fin >> n;
	for (int i = 1; i <= n; ++i)
		for (int j = 1; j <= n; ++j)
			fin >> nums[i][j];
	fin.close();
	//timer
	timespec_get(&sts, TIME_UTC);
	avx512_Gauss(n);
	timespec_get(&ets, TIME_UTC);
	print_result("avx512 gauss");

	for (int i = 1; i <= n; ++i) {
		for (int j = 1; j <= n; ++j)
			fout << nums[i][j] << " ";
		fout << endl;
	}
	fout.close();
	memset(nums, 0, sizeof nums);

}